#sig = sd(lc)
u = y - plogis(lc,location = -median(lc),scale = kc); n = length(u)
#u = y - plogis(1+x%*%c(b)); n = length(u)
#z = cbind(x,scale(x,scale = FALSE)^2)
z = x
n*energy::dcov(u,z)^2
}
fndcreg3(c(0.8,-0.8))
fndcreg3(c(1,-1))
(dcopt=optim(c(.5,-.5),fndcreg3))
kc = sqrt(3)/pi
print(dcopt$par/kc)
#print(dcopt$par)
} #end for
fndcreg3 = function(b){
lc = x%*%c(b); lc = lc/sd(lc)
#fn = function(t) {mean(y - plogis(t+))^2}
#robj=uniroot(f=fnt,c(-15,15)); th0=robj$root
#th0 = optimise(fn,c(-5,5))$minimum
#u = y - plogis(-median(lc)/kc+lc/kc); n = length(u)
#sig = sd(lc)
#u = y - plogis(lc,location = -median(lc),scale = kc); n = length(u)
u = y - plogis(lc,location = -mean(lc),scale = kc); n = length(u)
#u = y - plogis(1+x%*%c(b)); n = length(u)
#z = cbind(x,scale(x,scale = FALSE)^2)
z = x
n*energy::dcov(u,z)^2
}
fndcreg3(c(0.8,-0.8))
fndcreg3(c(1,-1))
(dcopt=optim(c(.5,-.5),fndcreg3))
kc = sqrt(3)/pi
print(dcopt$par/kc)
for (l in 1:20) {
set.seed(l)
x = matrix(rnorm(2*n),ncol = 2)
lc = 0 + x%*%c(1,-1) + rlogis(n)
y = as.numeric(I(lc>0))
fndcreg3 = function(b){
lc = x%*%c(b); lc = lc/sd(lc)
#fn = function(t) {mean(y - plogis(t+))^2}
#robj=uniroot(f=fnt,c(-15,15)); th0=robj$root
#th0 = optimise(fn,c(-5,5))$minimum
#u = y - plogis(-median(lc)/kc+lc/kc); n = length(u)
#sig = sd(lc)
#u = y - plogis(lc,location = -median(lc),scale = kc); n = length(u)
u = y - plogis(lc,location = -mean(lc),scale = kc); n = length(u)
#u = y - plogis(1+x%*%c(b)); n = length(u)
#z = cbind(x,scale(x,scale = FALSE)^2)
z = x
n*energy::dcov(u,z)^2
}
fndcreg3(c(0.8,-0.8))
fndcreg3(c(1,-1))
(dcopt=optim(c(.5,-.5),fndcreg3))
kc = sqrt(3)/pi
print(dcopt$par/kc)
#print(dcopt$par)
} #end for
fndcreg3 = function(b){
lc = x%*%c(b); lc = lc/sd(lc)
u = y - plogis(lc,location = -mean(lc)); n = length(u)
z = x
n*energy::dcov(u,z)^2
}
fndcreg3(c(0.8,-0.8))
fndcreg3(c(1,-1))
(dcopt=optim(c(.5,-.5),fndcreg3))
print(dcopt$par)
print(dcopt$par/kc)
kc
n=250
for (l in 1:20) {
set.seed(l)
x = matrix(rnorm(2*n),ncol = 2)
lc = 0 + x%*%c(1,-1) + rlogis(n)
y = as.numeric(I(lc>0))
fndcreg3 = function(b){
lc = x%*%c(b); lc = lc/sd(lc)
u = y - plogis(lc,location = -mean(lc),scale = kc); n = length(u)
z = x
n*energy::dcov(u,z)^2
}
fndcreg3(c(0.8,-0.8))
fndcreg3(c(1,-1))
(dcopt=optim(c(.5,-.5),fndcreg3))
kc = sqrt(3)/pi
print(dcopt$par/kc)
#print(dcopt$par)
} #end for
print(dcopt$par)
print(dcopt$par/kc)
dcopt$par
dcopt$par/kc
for (l in 1:20) {
set.seed(l)
x = matrix(rnorm(2*n),ncol = 2)
lc = 0 + x%*%c(1,-1) + rlogis(n)
y = as.numeric(I(lc>0))
fndcreg3 = function(b){
lc = x%*%c(b); lc = lc/sd(lc); kc = sqrt(3)/pi
u = y - plogis(lc,location = -mean(lc),scale = kc); n = length(u)
z = x
n*energy::dcov(u,z)^2
}
fndcreg3(c(0.8,-0.8))
fndcreg3(c(1,-1))
(dcopt=optim(c(.5,-.5),fndcreg3))
kc = sqrt(3)/pi
#print(dcopt$par/kc)
print(dcopt$par)
} #end for
for (l in 1:20) {
set.seed(l)
x = matrix(rnorm(2*n),ncol = 2)
lc = 0 + x%*%c(1,-1) + rlogis(n)
y = as.numeric(I(lc>0))
fndcreg3 = function(b){
lc = x%*%c(b); lc = lc/sd(lc); kc = sqrt(3)/pi
u = y - plogis(lc,location = -mean(lc),scale = kc); n = length(u)
z = x
n*energy::dcov(u,z)^2
}
fndcreg3(c(0.8,-0.8))
fndcreg3(c(1,-1))
(dcopt=optim(c(.5,-.5),fndcreg3))
kc = sqrt(3)/pi
print(dcopt$par/kc)
#print(dcopt$par)
} #end for
install.packages("expectreg v0.50")
install.packages("expectreg")
install.packages("erboost")
library(erboost)
vignette(erboost)
set.seed(12)
n = 200
er = rchisq(n,df=1)
x = mvrnorm(n=n,m=c(0,0),Sigma = matrix(c(1,0.25,0.25,1),ncol = 2))
set.seed(12)
n = 200
er = rchisq(n,df=1)
x = MASS::mvrnorm(n=n,m=c(0,0),Sigma = matrix(c(1,0.25,0.25,1),ncol = 2))
y = exp(-1 + x%*%c(-1,1)) + er
summary(cbind(x,y))
summary(data.frame(x,y))
thfun = function(th0,theta,y,x){
er=y - exp(th0+x%*%c(theta))
sum(er^2)
}
thfun = function(th0,theta,y,x){
elc = exp(x%*%c(theta))
log(sum(y*elc)/sum(elc^2))
}
thfun = function(theta,y,x){
elc = exp(x%*%c(theta))
log(sum(y*elc)/sum(elc^2))
}
thfun(theta=c(-1,1),y=y,x=x)
set.seed(12)
n = 2000
er = rchisq(n,df=1)
x = MASS::mvrnorm(n=n,m=c(0,0),Sigma = matrix(c(1,0.25,0.25,1),ncol = 2))
y = exp(-1 + x%*%c(-1,1)) + er
summary(data.frame(x,y))
thfun = function(theta,y,x){
elc = exp(x%*%c(theta))
log(sum(y*elc)/sum(elc^2))
}
thfun(theta=c(-1,1),y=y,x=x)
n = 2000
er = rchisq(n,df=1)
x = MASS::mvrnorm(n=n,m=c(0,0),Sigma = matrix(c(1,0.25,0.25,1),ncol = 2))
y = exp(-1 + x%*%c(-1,1)) + er
summary(data.frame(x,y))
thfun = function(theta,y,x){
elc = exp(x%*%c(theta))
log(sum(y*elc)/sum(elc^2))
}
thfun(theta=c(-1,1),y=y,x=x)
n = 200
er = rchisq(n,df=1)
x = MASS::mvrnorm(n=n,m=c(0,0),Sigma = matrix(c(1,0.25,0.25,1),ncol = 2))
y = exp(-1 + x%*%c(-1,1)) + er
summary(data.frame(x,y))
thfun = function(theta,y,x){
elc = exp(x%*%c(theta))
log(sum(y*elc)/sum(elc^2))
}
thfun(theta=c(-1,1),y=y,x=x)
n = 200
er = rchisq(n,df=1)
x = MASS::mvrnorm(n=n,m=c(0,0),Sigma = matrix(c(1,0.25,0.25,1),ncol = 2))
y = exp(-1 + x%*%c(-1,1)) + er
summary(data.frame(x,y))
thfun = function(theta,y,x){
elc = exp(x%*%c(theta))
log(sum(y*elc)/sum(elc^2))
}
thfun(theta=c(-1,1),y=y,x=x)
set.seed(12)
n = 200
er = rchisq(n,df=1)
er = rnorm(n)
x = MASS::mvrnorm(n=n,m=c(0,0),Sigma = matrix(c(1,0.25,0.25,1),ncol = 2))
y = exp(-1 + x%*%c(-1,1)) + er; y[y<0]=0
summary(data.frame(x,y))
thfun = function(theta,y,x){
elc = exp(x%*%c(theta))
log(sum(y*elc)/sum(elc^2))
}
thfun(theta=c(-1,1),y=y,x=x)
set.seed(12)
n = 200
er = rchisq(n,df=1)
#er = rnorm(n)
x = MASS::mvrnorm(n=n,m=c(0,0),Sigma = matrix(c(1,0.25,0.25,1),ncol = 2))
y = exp(-1 + x%*%c(-1,1)) + er; y[y<0]=0
summary(data.frame(x,y))
thfun = function(theta,y,x){
elc = exp(x%*%c(theta))
log(sum(y*elc)/sum(elc^2))
}
thfun(theta=c(-1,1),y=y,x=x)
set.seed(12)
n = 200
#er = rchisq(n,df=1)
er = rnorm(n)
x = MASS::mvrnorm(n=n,m=c(0,0),Sigma = matrix(c(1,0.25,0.25,1),ncol = 2))
y = exp(-1 + x%*%c(-1,1)) + er; y[y<0]=0
summary(data.frame(x,y))
thfun = function(theta,y,x){
elc = exp(x%*%c(theta))
log(sum(y*elc)/sum(elc^2))
}
thfun(theta=c(-1,1),y=y,x=x)
set.seed(12)
n = 200
#er = rchisq(n,df=1)
er = rnorm(n)
x = MASS::mvrnorm(n=n,m=c(0,0),Sigma = matrix(c(1,0.25,0.25,1),ncol = 2))
y = exp(-1 + x%*%c(-1,1)) + er; y[y<0]=0
summary(data.frame(x,y))
thfun = function(theta,y,x){
elc = exp(x%*%c(theta))
log(sum(y*elc)/sum(elc^2))
}
thfun(theta=c(-1,1),y=y,x=x)
thfun = function(th0,theta,y,x){
er =y- exp(th0+x%*%c(theta))
er = er-mean(er)
sum(er^2)
}
?optimise
thfun2 = function(th0,theta,y,x){
er =y- exp(th0+x%*%c(theta))
er = er-mean(er)
sum(er^2)
}
optimise(thfun2,c(-5,5),theta=c(-1,1),y=y,x=x)
thfun3 = function(th0,theta,y,x){
er =y- exp(th0+x%*%c(theta))
#er = er-mean(er)
sum(er^2)
}
optimise(thfun3,c(-5,5),theta=c(-1,1),y=y,x=x)
set.seed(12)
n = 200
#er = rchisq(n,df=1)
er = rnorm(n)
x = MASS::mvrnorm(n=n,m=c(0,0),Sigma = matrix(c(1,0.25,0.25,1),ncol = 2))
y = exp(-1 + x%*%c(-1,1)) + er; y[y<0]=0
summary(data.frame(x,y))
thfun = function(theta,y,x){
elc = exp(x%*%c(theta))
log(sum(y*elc)/sum(elc^2))
}
thfun(theta=c(-1,1),y=y,x=x)
thfun2 = function(th0,theta,y,x){
er =y- exp(th0+x%*%c(theta))
er = er-mean(er)
sum(er^2)
}
optimise(thfun2,c(-5,5),theta=c(-1,1),y=y,x=x)
thfun2 = function(th0,theta,y,x){
er =y- exp(th0+x%*%c(theta))
#er = er-mean(er)
sum(er^2)
}
optimise(thfun2,c(-5,5),theta=c(-1,1),y=y,x=x)
thfun3 = function(th0,theta,y,x){
er =y- exp(th0+x%*%c(theta))
er = er-mean(er)
sum(er^2)
}
optimise(thfun3,c(-5,5),theta=c(-1,1),y=y,x=x)
set.seed(12)
n = 200
er = rchisq(n,df=1)
#er = rnorm(n)
x = MASS::mvrnorm(n=n,m=c(0,0),Sigma = matrix(c(1,0.25,0.25,1),ncol = 2))
y = exp(-1 + x%*%c(-1,1)) + er; y[y<0]=0
summary(data.frame(x,y))
thfun = function(theta,y,x){
elc = exp(x%*%c(theta))
log(sum(y*elc)/sum(elc^2))
}
thfun(theta=c(-1,1),y=y,x=x)
thfun2 = function(th0,theta,y,x){
er =y- exp(th0+x%*%c(theta))
#er = er-mean(er)
sum(er^2)
}
optimise(thfun2,c(-5,5),theta=c(-1,1),y=y,x=x)
thfun3 = function(th0,theta,y,x){
er =y- exp(th0+x%*%c(theta))
er = er-mean(er)
sum(er^2)
}
optimise(thfun3,c(-5,5),theta=c(-1,1),y=y,x=x)
R = 2000; bmat = matrix(NA,R,2)
set.seed(12)
n = 200
for (l in 1:R) {
er = rchisq(n,df=2)
#er = rnorm(n)
x = MASS::mvrnorm(n=n,m=c(0,0),Sigma = matrix(c(1,0.25,0.25,1),ncol = 2))
y = exp(-1 + x%*%c(-1,1)) + er; y[y<0]=0
summary(data.frame(x,y))
thfun = function(theta,y,x){
elc = exp(x%*%c(theta))
log(sum(y*elc)/sum(elc^2))
}
bmat[l,1]=thfun(theta=c(-1,1),y=y,x=x)
# thfun2 = function(th0,theta,y,x){
#   er =y- exp(th0+x%*%c(theta))
#   #er = er-mean(er)
#   sum(er^2)
# }
# optimise(thfun2,c(-5,5),theta=c(-1,1),y=y,x=x)
thfun3 = function(th0,theta,y,x){
er =y- exp(th0+x%*%c(theta))
er = er-mean(er)
sum(er^2)
}
bmat[l,2]=optimise(thfun3,c(-5,5),theta=c(-1,1),y=y,x=x)$minimum
}
summary(bmat)
summary(bmat-1)
summary(1-bmat)
summary(bmat)
summary(bmat+1)
n
set.seed(12)
er = rlogis(n)
x = MASS::mvrnorm(n=n,m=c(0,0),Sigma = matrix(c(1,0.25,0.25,1),ncol = 2))
y = as.numeric(I((-1 + x%*%c(-1,1) + er) >0))
summary(data.frame(x,y))
thfunlg = function(th0,theta,y,x){
er =y- plogis(th0+x%*%c(theta))
er = er-mean(er)
sum(er^2)
}
set.seed(12)
er = rlogis(n)
x = MASS::mvrnorm(n=n,m=c(0,1),Sigma = matrix(c(1,0.25,0.25,1),ncol = 2))
y = as.numeric(I((-1 + x%*%c(-1,1) + er) >0))
summary(data.frame(x,y))
glm(y~x,family = "binomial")
binomial
(glgobj=glm(y~x,family = "binomial"))
thfunlg = function(th0,theta,y,x){
er =y- plogis(th0+x%*%c(theta))
er = er-mean(er)
sum(er^2)
}
(blg=optimise(thfunlg,c(-5,5),theta=c(-1,1),y=y,x=x)$minimum)
(blg2=optimise(thfunlg,c(-5,5),theta=glgobj$coefficients[-1],y=y,x=x)$minimum)
thfunlg = function(th0,theta,y,x){
er =y- plogis(th0+x%*%c(theta))
#er = er-mean(er)
sum(er^2)
}
(blg1=optimise(thfunlg,c(-5,5),theta=c(-1,1),y=y,x=x)$minimum)
(blg2=optimise(thfunlg,c(-5,5),theta=glgobj$coefficients[-1],y=y,x=x)$minimum)
2.39/40.945
round(2.39/40.945,2)
round(100*c(0.6626016, 0.6618257, 0.9660000, 0.8840000),2)
eps = rnorm(1000)
Z = abs(rnorm(1000)+eps)
cor(eps,Z)
cor.test(eps,Z)
cor.test(eps,Z,method = "kendall")
energy::dcov.test(eps,Z,R=999)
energy::dcov.test(sign(eps),Z,R=999)
X=rnorm(1000);Z = abs(X+eps)
cor(eps,Z)
cor.test(eps,Z)
energy::dcov.test(eps,Z,R=999)
energy::dcov.test(sign(eps),Z,R=999)
energy::dcov.test(X,Z,R=999)
curve(abs,-2,2)
fn = function(x) x^(2/3)
curve(fn,-2,2,add = TRUE)
(-2)^(2/3)
fn = function(x) (x^2)^(1/3)
curve(abs,-2,2)
curve(fn,-2,2,add = TRUE,lty=2)
fn = function(x) (x^2)^(1/1.9)
curve(fn,-2,2,add = TRUE,lty=2)
curve(abs,-2,2);curve(fn,-2,2,add = TRUE,lty=2)
fn = function(x) (x^2)^(1/1.5)
curve(abs,-2,2);curve(fn,-2,2,add = TRUE,lty=2)
x=rnorm(1000)
cor.test(x,abs(x))
cor.test(x,fn(x))
fn = function(x) (x^2)^(1/1.75)
curve(abs,-2,2);curve(fn,-2,2,add = TRUE,lty=2)
cor.test(x,fn(x))
cor.test(x,abs(x))
er = (rchisq(1000,1)-1)/sqrt(2)
summary(er)
plot(density(er))
cor.test(er,abs(er))
plot(density(er))
cor.test(er,er^4)
summary(er)
er=er-mean(er)
summary(er)
cor.test(er,er^4)
rd = rnorm(1000,mean = 0,sd=10)
summary(rd)
rd = rnorm(1000,mean = 0,sd=3); summary(rd)
?dist
library(MASS)
library(mvtnorm)
library(GJRM)
library(copula)
library(quantreg)
library(MASS)
library(mvtnorm)
library(GJRM)
library(doParallel)
library(doRNG)
cl <- makeCluster(2, type = "FORK")
registerDoParallel(cl)
registerDoRNG(seed = 123)
foreach(i=1:4, .combine = 'c') %dopar% {print(i);rnorm(1, mean = 0, sd = 1)}
stopImplicitCluster()
library(doParallel)
library(doRNG)
cl <- makeCluster(2, type = "FORK")
registerDoParallel(cl)
registerDoRNG(seed = 123)
foreach(i=1:4, .combine = 'c') %dopar% {message("i = ",i)
rnorm(1, mean = 0, sd = 1)}
stopImplicitCluster()
?mosaic::set.rseed
library(EDMeasure)
n = 200
set.seed(n)
x=rchisq(n,df=1)
x = scale(x,scale = FALSE)
Tx = abs(x)
summary(cbind(x,Tx))
cor.test(x,Tx)
xTx.n=cbind(x,Tx)%*%solve(expm::sqrtm(cov(cbind(x,Tx))))
cor(xTx.n)
cor.test(xTx.n[,1],xTx.n[,2])
energy::dcorT.test(xTx.n[,1],xTx.n[,2]) #still dependent but uncorrelat
EDMeasure::mdc(xTx.n[,1],xTx.n[,2])
fntrans<- function(a,X){
TX = abs(X-a) #transformed X
X = scale(X,scale = FALSE)
xTx.n = cbind(X,TX)
xTx.n = xTx.n%*%solve(expm::sqrtm(cov(xTx.n)))
EDMeasure::mdc(xTx.n[,1],xTx.n[,2])
}
set.seed(n)
x=rchisq(n,df=1)
summary(x)
fntrans(a=mean(x),X=x)
fntrans(a=median(x),X=x)
fn = function(a)fntrans(a=a,X=x)
fn = Vectorize(fn)
curve(fn,0.2,mean(x))
curve(fn,from=0.2,to=mean(x))
(opta=(optimise(fn,c(0.2,mean(x)))))
TX = abs(x-opta$minimum) #transformed X
X = scale(x,scale = FALSE) #demeaned X
xTx.n = cbind(X,TX)
xTx.n = xTx.n%*%solve(expm::sqrtm(cov(xTx.n))) #remove correlation
EDMeasure::cmdm_test(xTx.n[,1],xTx.n[,2],xTx.n[,2]) #fail to reject!
energy::dcorT.test(xTx.n[,1],xTx.n[,2])
#===================================================================================>
#This file saves the Stata dataset to txt format
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) #set working directory path to source file
library(foreign)
#===================================================================================>
list.files()
#===================================================================================>
#This file saves the Stata dataset to txt format
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) #set working directory path to source file
library(foreign)
#===================================================================================>
list.files()
hornung_data_textiles<- read.dta("hornung-data-textiles.dta")
dim(hornung_data_textiles)
hornung_data_textiles[1,]
# save file into txt format
write.table(hornung_data_textiles, "hornung_data_textiles.txt", sep = " ",col.names = TRUE)
list.files()
